#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

          //////////////////
          // Input   ///////
          //////////////////

	  int Nf;
	  int NOa2;
	  double ainv;
          char dir[1024];

    
          ////// Parsing command line

          for(int i=1; i<argc; i++) {
            if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
            if(strncmp(argv[i],"-Nf",3)==0)  sscanf(argv[i],"-Nf=%d",&Nf);
            if(strncmp(argv[i],"-dir",4)==0)  sscanf(argv[i],"-dir=%s",dir);
            if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T);
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L);
            if(strncmp(argv[i],"-ainv",5)==0)  sscanf(argv[i],"-ainv=%lf",&ainv);
            if(strncmp(argv[i],"-NOa2",5)==0)  sscanf(argv[i],"-NOa2=%d",&NOa2);
          }

   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
  read_momenta_data(momentafile);
  Equivalent_Momenta(SinP2);
  double *nn=(double*)malloc(sizeof(double)*(Nmom));
  compute_nn(nn);

  // IO

   char file[1024];
   char output[1024];
   char Oa2print[1024];
   char dir_qav[1024];
   sprintf(dir_qav,"%s/QAV/",dir);

    //////////////////////////////////
    // Read Zq 
     /////////////////////////////

    double ***Zq=malloc(sizeof(double**)*(Nmom_eq));
    for (int imom=0; imom<Nmom_eq; imom++) Zq[imom]=malloc(sizeof(double*)*(NOa2));
    for (int imom=0; imom<Nmom_eq; imom++) for(int iOa2=0; iOa2<NOa2; iOa2++) Zq[imom][iOa2]=malloc(sizeof(double)*(Nclust+1));

    O55 ***ZpPC=malloc(sizeof(O55**)*(Nmom_eq));
    for (int imom=0; imom<Nmom_eq; imom++) ZpPC[imom]=malloc(sizeof(O55*)*(NOa2));
    for (int imom=0; imom<Nmom_eq; imom++) for(int iOa2=0; iOa2<NOa2; iOa2++) ZpPC[imom][iOa2]=malloc(sizeof(O55)*(Nclust+1));

/*  O55 ***DeltapPC=malloc(sizeof(O55**)*(Nmom_eq));
    for (int imom=0; imom<Nmom_eq; imom++) DeltapPC[imom]=malloc(sizeof(O55*)*(NOa2));
    for (int imom=0; imom<Nmom_eq; imom++) for(int iOa2=0; iOa2<NOa2; iOa2++) DeltapPC[imom][iOa2]=malloc(sizeof(O55)*(Nclust+1));
*/
    double *aux_clust=malloc(sizeof(double)*(Nclust+1));
 
    for (int iOa2=0; iOa2<NOa2; iOa2++)
    {
    whichOa2(iOa2,Oa2print);
    sprintf(file,"Zq_%s",Oa2print);
    for (int imom=0; imom<Nmom_eq; imom++) {
            read_input_onemomenta(dir_qav,file,aux_clust,imom);
            for (int iclust=0; iclust<Nclust+1; iclust++) Zq[imom][iOa2][iclust]=aux_clust[iclust];
     }

    for (int imom=0; imom<Nmom_eq; imom++) { 
        for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
          {
                sprintf(file,"ZpPC%d%d_%s",iop1+1,iop2+1,Oa2print);
                read_input_onemomenta(dir_qav,file,aux_clust,imom);
                for (int iclust=0; iclust<Nclust+1; iclust++) ZpPC[imom][iOa2][iclust][iop1][iop2]=aux_clust[iclust];
          }
    }
/*
    for (int imom=0; imom<Nmom_eq; imom++) {
        for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
          { 
                sprintf(file,"DeltapPC%d%d_%s",iop1+1,iop2+1,Oa2print);
                read_input_onemomenta(dir_qav,file,aux_clust,imom);
                for (int iclust=0; iclust<Nclust+1; iclust++) DeltapPC[imom][iOa2][iclust][iop1][iop2]=aux_clust[iclust];
          }
    }
*/
   }

   ///////////////////////////////////////
   ////// COMPUTE ZRGI
   ///////////////////////////////////////

    double ***Zq_RGI=malloc(sizeof(double**)*(Nmom_eq));
    for (int imom=0; imom<Nmom_eq; imom++) Zq_RGI[imom]=malloc(sizeof(double*)*(NOa2));
    for (int imom=0; imom<Nmom_eq; imom++) for(int iOa2=0; iOa2<NOa2; iOa2++) Zq_RGI[imom][iOa2]=malloc(sizeof(double)*(Nclust+1));

    O55 ***ZpPC_RGI=malloc(sizeof(O55**)*(Nmom_eq));
    for (int imom=0; imom<Nmom_eq; imom++) ZpPC_RGI[imom]=malloc(sizeof(O55*)*(NOa2));
    for (int imom=0; imom<Nmom_eq; imom++) for(int iOa2=0; iOa2<NOa2; iOa2++) ZpPC_RGI[imom][iOa2]=malloc(sizeof(O55)*(Nclust+1));

    for (int iOa2=0; iOa2<NOa2; iOa2++)
	for(int iclust=0; iclust<Nclust+1; iclust++)
         for (int imom=0; imom<Nmom_eq; imom++)
         {
		Zq_RGI[imom][iOa2][iclust]=evolve_Zq_to_RGI(Zq[imom][iOa2][iclust],SinP2_eq[imom],Nf,ainv);
		evolve_Z4F_to_RGI(ZpPC_RGI[imom][iOa2][iclust],ZpPC[imom][iOa2][iclust],SinP2_eq[imom],Nf,ainv);
	 }



        /////////////////
        // Write output

     double **aux_allmomenta_real=malloc(sizeof(double*)*(Nmom_eq));
     for (int imom=0; imom<Nmom_eq; imom++) aux_allmomenta_real[imom]=malloc(sizeof(double)*(Nclust+1));

   for (int iOa2=0; iOa2<NOa2; iOa2++)
   {
         whichOa2(iOa2,Oa2print);
         sprintf(output,"%s/Z_ALLP/Zq_RGI_%s.out",dir_qav,Oa2print);
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=Zq_RGI[imom][iOa2][iclust];
         write_output_allmomenta_eq_real(output,aux_allmomenta_real,SinP2_eq);

         sprintf(file,"Zq_RGI_%s",Oa2print);
         write_output_momentabymomenta_eq_real(dir_qav,file,aux_allmomenta_real,SinP2_eq);
   }

   for (int iOa2=0; iOa2<NOa2; iOa2++)
   {
        whichOa2(iOa2,Oa2print);
        for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
        {
         sprintf(output,"%s/Z_ALLP/ZpPC%d%d_RGI_%s.out",dir_qav,iop1+1,iop2+1,Oa2print);
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=ZpPC_RGI[imom][iOa2][iclust][iop1][iop2];
         write_output_allmomenta_eq_real(output,aux_allmomenta_real,SinP2_eq);

         sprintf(file,"ZpPC%d%d_RGI_%s",iop1+1,iop2+1,Oa2print);
         write_output_momentabymomenta_eq_real(dir_qav,file,aux_allmomenta_real,SinP2_eq);
        }
   }
/*   for (int iOa2=0; iOa2<NOa2; iOa2++)
   {
        whichOa2(iOa2,Oa2print);
        for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
        {
         sprintf(output,"%s/Z_ALLP/DeltapPC%d%d_%s.out",dir_qav,iop1+1,iop2+1,Oa2print);
         for (int imom=0; imom<Nmom_eq; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta_real[imom][iclust]=DeltapPC[imom][iOa2][iclust][iop1][iop2];
         write_output_allmomenta_eq_real(output,aux_allmomenta_real,SinP2_eq);

         sprintf(file,"DeltapPC%d%d_%s",iop1+1,iop2+1,Oa2print);
         write_output_momentabymomenta_eq_real(dir_qav,file,aux_allmomenta_real,SinP2_eq);
        }
   }

*/
    free(Zq); free(ZpPC); free(Zq_RGI); free(ZpPC_RGI);
        ////////////////

    return 0;

}

